3  ANOVA_oneway

Author

Brayan Cubides

4 Introducción

Este documento presenta un análisis estadístico en dos partes. La primera parte explora la relación entre dos variables categóricas (hábito de fumar de la madre y bajo peso del bebé al nacer) utilizando tablas de contingencia y la prueba \(\chi^2\) de independencia. La segunda parte introduce el Análisis de Varianza (ANOVA) de una vía para comparar las medias de una variable cuantitativa (peso de plantas) a través de diferentes grupos de tratamiento.

4.0.1 Librerías Requeridas

Se cargan los paquetes necesarios para el análisis.

library(readxl)
library(dplyr)
library(ggplot2)
library(RcmdrMisc) # Para análisis de tablas de contingencia
library(gplots)    # Para plotmeans()
library(ggpubr)
library(multcomp)  # Para comparaciones múltiples
library(car)       # Para Levene's Test

4.0.2 Carga de Datos

pto9 <- read_excel("dataset_trabajo1.xlsx")

5 Parte 1: Test Chi-Cuadrado de Independencia

Se investiga si existe una relación estadísticamente significativa entre el hábito de fumar de la madre (smoker) y la probabilidad de que el bebé nazca con bajo peso (lowbwt).

5.1 Tablas de Contingencia

Estas tablas resumen las frecuencias de las combinaciones de las dos variables.

  • Frecuencias Absolutas: Conteos directos.
  • Perfiles Fila/Columna: Proporciones condicionales que ayudan a visualizar la relación.
# Frecuencias Absolutas con totales
addmargins(table(pto9$smoker, pto9$lowbwt))
     
       0  1 Sum
  0   19  1  20
  1   17  5  22
  Sum 36  6  42
# Análisis más detallado con xtabs()
Table <- xtabs(~ pto9$smoker + pto9$lowbwt)
addmargins(Table)      # Frecuencias Absolutas
           pto9$lowbwt
pto9$smoker  0  1 Sum
        0   19  1  20
        1   17  5  22
        Sum 36  6  42
totPercents(Table)     # Frecuencias Relativas (porcentaje del total)
         0    1 Total
0     45.2  2.4  47.6
1     40.5 11.9  52.4
Total 85.7 14.3 100.0
rowPercents(Table)     # Perfiles Fila
           pto9$lowbwt
pto9$smoker    0    1 Total Count
          0 95.0  5.0   100    20
          1 77.3 22.7   100    22
colPercents(Table)     # Perfiles Columna
           pto9$lowbwt
pto9$smoker     0     1
      0      52.8  16.7
      1      47.2  83.3
      Total 100.0 100.0
      Count  36.0   6.0

5.2 Gráfico de Frecuencias

Un gráfico de barras apiladas permite visualizar las frecuencias y proporciones.

par(mfrow = c(1, 2))
# Gráfico de Frecuencias Absolutas
ggplot(pto9, aes(x = as.factor(smoker), fill = as.factor(lowbwt))) +
  geom_bar() +
  labs(title = "Frecuencias Absolutas", x = "Madre Fumadora", y = "Conteo", fill = "Bajo Peso") +
  theme_minimal()

# Gráfico de Perfiles (Proporciones)
ggplot(pto9, aes(x = as.factor(smoker), fill = as.factor(lowbwt))) +
  geom_bar(position = "fill") +
  labs(title = "Perfiles (Proporciones)", x = "Madre Fumadora", y = "Proporción", fill = "Bajo Peso") +
  theme_minimal()

5.3 Test Chi-Cuadrado (\(\chi^2\)) de Independencia

Este test evalúa formalmente la hipótesis de independencia entre las dos variables.

  • \(H_0\): Las variables son independientes.
  • \(H_1\): Las variables no son independientes (están asociadas).
  • Estadístico: \(\chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}}\)
# Test estándar
chi <- chisq.test(as.factor(pto9$smoker), as.factor(pto9$lowbwt))
Warning in chisq.test(as.factor(pto9$smoker), as.factor(pto9$lowbwt)):
Chi-squared approximation may be incorrect
chi

    Pearson's Chi-squared test with Yates' continuity correction

data:  as.factor(pto9$smoker) and as.factor(pto9$lowbwt)
X-squared = 1.4358, df = 1, p-value = 0.2308
# Test con simulación de Monte Carlo (recomendado cuando hay frecuencias esperadas bajas)
chisq.test(as.factor(pto9$smoker), as.factor(pto9$lowbwt), simulate.p.value = TRUE, B = 2000)

    Pearson's Chi-squared test with simulated p-value (based on 2000
    replicates)

data:  as.factor(pto9$smoker) and as.factor(pto9$lowbwt)
X-squared = 2.6886, df = NA, p-value = 0.1859

Conclusión: Un p-valor > 0.05 sugiere que no hay evidencia suficiente para rechazar la hipótesis nula de independencia.

5.4 V de Cramer

Mide la fuerza de la asociación (0 = sin asociación, 1 = asociación perfecta).

V <- sqrt(chi$statistic / (length(pto9$smoker) * min(2 - 1, 2 - 1)))
V
X-squared 
0.1848935 

6 Parte 2: ANOVA de una Vía

Se utiliza el Análisis de Varianza (ANOVA) para comparar las medias de una variable cuantitativa (weight) entre tres o más grupos definidos por una variable categórica (group).

  • \(H_0\): Las medias de todos los grupos son iguales (\(\mu_{ctrl} = \mu_{trt1} = \mu_{trt2}\)).
  • \(H_1\): Al menos una de las medias es diferente.

6.1 1. Lectura de Datos y Exploración

Se utiliza el conjunto de datos PlantGrowth.

data <- PlantGrowth
data$group <- factor(data$group, levels = c("ctrl", "trt1", "trt2"))
head(data)
  weight group
1   4.17  ctrl
2   5.58  ctrl
3   5.18  ctrl
4   6.11  ctrl
5   4.50  ctrl
6   4.61  ctrl

6.2 2. Resúmenes Estadísticos por Grupo

group_by(data, group) %>%
  summarise(
    count = n(),
    mean = mean(weight, na.rm = TRUE),
    sd = sd(weight, na.rm = TRUE)
  )
# A tibble: 3 × 4
  group count  mean    sd
  <fct> <int> <dbl> <dbl>
1 ctrl     10  5.03 0.583
2 trt1     10  4.66 0.794
3 trt2     10  5.53 0.443

6.3 3. Diagnóstico Gráfico

Se visualizan los datos para evaluar si los supuestos del ANOVA (normalidad, homocedasticidad) parecen razonables.

# Boxplots
ggboxplot(data, x = "group", y = "weight", 
          color = "group", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
          order = c("ctrl", "trt1", "trt2"),
          ylab = "Peso", xlab = "Tratamiento")

# Gráfico de Medias con Intervalos de Confianza
ggline(data, x = "group", y = "weight", 
       add = c("mean_se", "jitter"), 
       order = c("ctrl", "trt1", "trt2"),
       ylab = "Peso", xlab = "Tratamiento")

6.4 4. Aplicación del Modelo ANOVA

Se ajusta el modelo ANOVA para obtener la tabla de análisis de varianza y la prueba F.

# Estimación de los efectos (tau)
hat_mu <- mean(data$weight)
tau_ctrl <- mean(data$weight[data$group == "ctrl"]) - hat_mu
tau_tr1 <- mean(data$weight[data$group == "trt1"]) - hat_mu
tau_tr2 <- mean(data$weight[data$group == "trt2"]) - hat_mu
c(tau_ctrl, tau_tr1, tau_tr2)
[1] -0.041 -0.412  0.453
# Ajuste del modelo ANOVA
res.aov <- aov(weight ~ group, data = data)
summary(res.aov)
            Df Sum Sq Mean Sq F value Pr(>F)  
group        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusión: El p-valor (0.0159) es menor que 0.05, por lo que se rechaza la hipótesis nula. Se concluye que existe una diferencia estadísticamente significativa en el peso promedio entre al menos dos de los grupos de tratamiento.

6.5 5. Comparaciones por Pares (Post-Hoc)

Una vez que se rechaza la \(H_0\) global, se realizan pruebas post-hoc para identificar qué pares de grupos son significativamente diferentes.

# Método de Tukey HSD (Honest Significant Difference)
TukeyHSD(res.aov)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = weight ~ group, data = data)

$group
            diff        lwr       upr     p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1  0.865  0.1737839 1.5562161 0.0120064
# Pruebas t por pares con corrección de Bonferroni
pairwise.t.test(data$weight, data$group, p.adjust.method = "bonferroni")

    Pairwise comparisons using t tests with pooled SD 

data:  data$weight and data$group 

     ctrl  trt1 
trt1 0.583 -    
trt2 0.263 0.013

P value adjustment method: bonferroni 

6.6 6. Validación de Supuestos del Modelo ANOVA

Se verifican los supuestos del modelo sobre los residuales.

6.6.1 a. Homocedasticidad (Igualdad de Varianzas)

  • \(H_0\): Las varianzas de los grupos son iguales.
# Gráfico de diagnóstico
plot(res.aov, 1)

# Test de Levene (más robusto)
leveneTest(weight ~ group, data = data)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2  1.1192 0.3412
      27               
# Test de Bartlett (sensible a la no normalidad)
bartlett.test(weight ~ group, data = data)

    Bartlett test of homogeneity of variances

data:  weight by group
Bartlett's K-squared = 2.8786, df = 2, p-value = 0.2371

Conclusión: Los p-valores de ambos tests son grandes, por lo que no se rechaza la hipótesis de homocedasticidad.

6.6.2 b. Normalidad de los Residuales

  • \(H_0\): Los residuales siguen una distribución normal.
# Gráfico Q-Q Plot
plot(res.aov, 2)

# Test de Shapiro-Wilk
shapiro.test(res.aov$residuals)

    Shapiro-Wilk normality test

data:  res.aov$residuals
W = 0.96607, p-value = 0.4379

Conclusión: El p-valor es grande, por lo que no se rechaza el supuesto de normalidad.

6.7 7. Alternativa No Paramétrica: Test de Kruskal-Wallis

Si los supuestos del ANOVA no se cumplen (especialmente normalidad y homocedasticidad), se puede utilizar la prueba de Kruskal-Wallis, que es su análoga no paramétrica y trabaja con los rangos de los datos.

kruskal.test(weight ~ group, data = data)

    Kruskal-Wallis rank sum test

data:  weight by group
Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842